infile = "./Data_2.csv"

data <- read.csv(file=infile,head=TRUE,sep=",")
# heads are "Momentum" and "Position" and "Variance"


x = data$Position
p = data$Momentum
V = data$Variance

# Invert the data to transform a nonlinear fitting problem (Lorentzian) into a linear one (Parabola)

y = 1 / p

# Fit the data

lorentzian <- lm( y ~ x + I((x)^2) )

# Grab the coefficients and their std errors

alpha <- lorentzian$coefficients[1]
beta  <- lorentzian$coefficients[2]
delta <- lorentzian$coefficients[3]

print(c(alpha,beta,delta))

SE <- coef(summary(lorentzian))[,2]
# SE[1] is sigma_alpha, etc.

# Construct th physically meaningful coefficients

x_zero <- - beta / (2*delta)
print(x_zero)

gamma_squared <- alpha/delta - x_zero^2
g <- gamma_squared
print(g)
gamma <- sqrt(g)
print(gamma)

I <- 1/(gamma_squared*delta)
print(I)

# Construct the Prediction Interval

sigma_x_zero <- abs(x_zero)/2 * sqrt((SE[2]/beta)^2+(SE[3]/delta)^2)
temp <- (SE[1]/alpha)^2+(SE[3]/delta)^2+2*abs(x_zero)*sigma_x_zero
sigma_g <- gamma_squared * sqrt(temp)
sigma_I <- I^2 * sqrt((SE[3]/delta)^2+(sigma_g/g)^2)

print(sigma_x_zero)
print(temp)
print(sigma_g)
print(sigma_I)

fitted <- function(x) { I*gamma_squared / ((x - x_zero)^2 + gamma_squared)}

sigma_fitted <- function(x) {
    a <- ( (2*abs(x_zero)*sigma_x_zero)^2 + (sigma_g)^2 )
    b <- ((x - x_zero)^2 + gamma_squared)^2
    c <- (SE[3]/delta)^2
    fitted(x)*sqrt(a/b+c)
}

print(sigma_fitted(x)[3:17])

upper = fitted(x)+2*sigma_fitted(x)
lower = fitted(x)-2*sigma_fitted(x)

# Plot it!

plot(x,p,pch=".",xlab="Position (Bohr)", ylab="Conditional Momentum (Atomic Units)")
points(x,upper, pch=".", col="blue",cex=0.8)
points(x,lower, pch=".", col="blue",cex=0.8)
curve(fitted, -0.7,-0.05, add=TRUE,col="dark violet")